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ABSTRACT 

We study the properties of cosmological shock waves identified in high- 
resolution, N-body/hydrodynamic simulations of a ACDM universe and their 
role on thermalization of gas and acceleration of nonthermal, cosmic ray (CR) 
particles. External shocks form around sheets, filaments and knots of mass distri- 
bution when the gas in void regions accretes onto them. Within those nonlinear 
structures, internal shocks are produced by infall of previously shocked gas to 
filaments and knots, and during subclump mergers, as well as by chaotic fiow 
motions. Due to the low temperature of the accreting gas, the Mach number 
of external shocks is high, extending up to M ~ 100 or higher. In contrast, 
internal shocks have mostly low Mach numbers. For all shocks of M > 1.5 the 
mean distance between shock surfaces over the entire computed volume is ~ Ah~^ 
Mpc at present, or ~ lh~^ Mpc for internal shocks within nonlinear structures. 
Identified external shocks are more extensive, with their surface area ~ 2 times 
larger than that of identified internal shocks at present. However, especially be- 
cause of higher preshock densities, but also due to higher shock speeds, internal 
shocks dissipate more energy. Hence, the internal shocks are mainly responsi- 
ble for gas thermalization as well as CR acceleration. In fact, internal shocks 
with 2 < M < 4 contribute ~ 1/2 of the total dissipation. Using a nonlinear 
diffusive shock acceleration model for CR protons, we estimate the ratio of CR 
energy to gas thermal energy dissipated at cosmological shock waves to be ~ 1/2 
through the history of the universe. Our result supports scenarios in which the 
intracluster medium contains energetically significant populations of CRs. 

Subject headings: large-scale structure of universe - methods:numerical - shock 
waves 
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1. Introduction 

According to cosmological N-body/hydrodynamic simulations, intergalactic shock waves 
develop as a consequence of the large scale structure formation of the universe. Infall of bary- 
onic gas toward sheets, filaments and knots, as well as supersonic flows associated with hierar- 
chical clustering, induce shocks (see, e.g., Quihs et al. 1998; Miniati et al. 2000; Miniati 2002; 
Gabici & Blasi 2003). Those cosmological shock waves, like most astrophysical shocks, are 
"coUisionless" features mediated by collective, electromagnetic viscosities. Such viscosities 
rely on irregular magnetic field components, i.e., MHD wave turbulence that is self-excited 
by the streaming suprathermal particles produced during shock formation (see, e.g., Wentzel 
1974; Kennel et al. 1985; Quest 1988). 

The existence and character of these shocks is important for several reasons. Through 
dissipation the cosmological shock waves convert part of the gravitational energy associated 
with structure formation into heat and, thus, govern the gas thermal history in the universe. 
Thermal energy is then radiated away, manifesting the large scale structure as well as the 
dynamics that created it (see, e.g., Cen & Ostriker 1999; Dave et al. 2001). At the same 
time, due to incomplete plasma thermalization at coUisionless shocks, a sizeable portion of 
the shock energy can be converted into cosmic ray (CR) energy (mostly ionic) via diffusive 
shock acceleration (DSA) (for reviews, see, e.g., Blandford & Eichler 1987; Malkov & Drury 
2001). This nonthermal population represents a small fraction of the ion flux through a 
coUisionless shock that leaks back upstream to become subject to DSA; that is, to be "in- 
jected" into the nonthermal population. Numerous nonlinear studies of DSA have shown 
that substantial fractions of the energy flux through strong shocks can be captured by the 
nonthermal populations (e.g., Berezhko et al. 1995; Ellison et al. 1996; Malkov 1999; Kang 
et al. 2002). Extensive nonlinear simulations by some of us incorporating a plasma-physics- 
based "thermal leakage" injection model into combined gas dynamic/CR diffusion-convection 
simulations found that strong shocks transfer up to ~ 1/2 of the initial shock kinetic energy 
to CRs by this process (Kang et al. 2002; Kang & Jones 2002). 

There is clear evidence that one or more processes energize signiflcant nonthermal parti- 
cle populations in and around cosmic structures. A number of clusters have been found with 
diffuse synchrotron radio halos or/and radio relic sources, indicating the existence of rela- 
tivistic electron populations in intracluster medium (ICM) (see, e.g., Giovannini & Feretti 
2000; Feretti 2003). In addition, some clusters have been reported to possess excess EUV 
and/or hard X-ray radiation compared to that expected from the hot, thermal X-ray emitting 
ICM, most likely produced by inverse Compton scattering of cosmic microwave background 
radiation (CMBR) photons by CR electrons (see, e.g.. Lieu et al. 1996; Fusco-Femiano et 
al. 1999). Also it has been suggested that a fraction of the diffuse 7-ray background radia- 
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tion could originate from the same process (Loeb & Waxmann 2000; Miniati 2002; Scharf & 
Mukherjee 2002). If some of those CR electrons have been energized at cosmological shock 
waves, the same process should have produced greater CR proton populations. Hence, al- 
though CR protons in the ICM have yet to be confirmed by the observation of 7-ray photons 
produced by inelastic collisions between CR and thermal-gas protons (sec, e.g., Miniati et 
al. 2001b; Reimer et al. 2003), there may very well exist CR proton populations there whose 
pressure is comparable to the gas thermal pressure (Enfihn et al. 1997; Colafrancesco & Blasi 
1998; Lieu et al. 1999). 

The properties of cosmological shock waves in the large scale structure of the universe 
were analyzed quantitatively by Miniati et al. (2000) using numerical simulations with 270^ 
grid zones in a cubic comoving region of size 85/i~^ Mpc for a SCDM model and a ACDM 
model. They identified accretion shocks, merger shocks and internal fiow shocks which were 
formed by infall and hierarchical clustering, and showed that the topology of these shocks is 
very complex. They found merger and internal fiow shocks distributed over Mach numbers 
from 3 to 10 with a peak at M ~ 5 and accretion shock Mach numbers ranging between 10 
and a few xlO^. In a recent study based on a higher resolution simulation with 512^ grid 
zones in a cubic comoving region of size 50/i-^ Mpc for a ACDM model, Miniati (2002) showed 
that most dissipation involves shocks of modest strength, with 4 < M < 10 accounting for 
~ 45% of total shock heating. On the other hand, through a semi-analytic study Gabici & 
Blasi (2003) found a Mach number distribution of merger-related shocks during large scale 
structure formation with a peak at much lower Mach number (M < 1.5). 

In this paper, we critically re-examine the properties of cosmological shock waves with 
a new set of cosmic structure formation simulations. We quantify the characteristics of cos- 
mological shock waves and estimate the dissipation, gas thermalization and CR acceleration 
at those shocks. The capture of shocks in hydrodynamic simulations and the identifica- 
tion of shocks in such simulation data are affected by numerical details including resolution. 
Therefore, to validate our findings we estimate the errors in our measured quantities through 
consistency checks and resolution convergence tests. This study adds valuable insights into 
the thermal history and nature of gas in the universe, as well as nonthermal activities in the 
ICM. In §2, simulations are detailed along with shock identification. The main results of 
shock characteristics and shock dissipation are described in §3 and §4, respectively, followed 
by a summary in §5. 
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2. 



Numerics 



2.1. 



Simulations 



The cold dark matter cosmology with a cosmological constant (ACDM) was employed 
with the following parameters: ^bm — 0.043, 0,dm — 0.227, and Qa — 0.73, h = ifo/(100 
km/s/Mpc) = 0.7, and as — 0.8. The above values are consistent with those fitted with 
the recent WMAP data (sec, e.g., Bennett et al. 2003). A cubic region of comoving size 
lOOhr^ Mpc was simulated inside a computational box with 1024'^, 512^, 256"^, 128'^ and 
64'^ grid zones for gas and gravity and 512'\ 256'^, 128^, 64"^ and 32"^ particles for dark 
matter. It allows a uniform spatial resolution of Al = 97. 7h^^ kpc — 1.57h~^ Mpc. The 
simulations were performed using a PM/Eulcrian hydrodynamic cosmology code. The code 
is described in Ryu et al. (1993). But the version of the code used includes several updates. 
For instance, it now adopts the MC (monotonized central difference) hmiter, instead of the 
original minmod limiter. The update to the MC limiter was intended to enhance the density 
resolution and to capture shocks more sharply (see, e.g., LeVeque 1997, for details). 

We did not include in our simulations several physical processes such as radiative cooling, 
galaxy/star formation, feedback from galaxies and stars, that can play significant roles in 
determining conditions within cluster cores, nor reionization of the intergalactic medium that 
effectively sets a temperature floor to the IGM. Our primary goal is to study cosmological 
shocks which are mostly outside cluster core regions. We established a temperature floor as 
a part of our analysis procedure. The conclusions drawn in this study, hence, should not be 
significantly weakened by the exclusion of those additional physical processes. 



While shocks are automatically detected during the simulations by the Riemann solver 
within the hydrodynamics routine, there are additional steps necessary to identify and char- 
acterize shocks for analysis. Wc have done this as a post-processing step using the simulation 
data at selected epochs. Ideally, explicitly three-dimensional flow motions should be con- 
sidered in identifying shocks in simulation data. However, to simplify the analysis we used 
a one-dimensional procedure applied successively in all three primary directions. A zone 
was tagged as a "shock zone" currently experiencing shock dissipation whenever these three 
criteria are met: 1) AT ■ As > 0, i.e., the gradients of gas temperature and entropy have 
the same sign, 2) V ■ v < 0, i.e., the local flow is converging (where V • v is the divergence 
of three-dimensional velocity field), 3) | A logT| > 0.11, where in each case we define central 
differences according to the scheme, AQ — Qi+i — Qi-i- The third condition corresponds 



2.2. 



Shock Identification 
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to the temperature jump of a Mach 1.3 shock. Typically, a shock is represented by a jump 
spread over 2 — 3 tagged zones. Hence, we identified a "shock center" within the numerical 
shock where V • v is minimum and labeled this center as part of a shock surface. 

The Mach number of shock centers, M, was calculated from the temperature jump 
across shocks, which is given by 

T2 ^ (5M^-l)(M^ + 3) 

Ti 16M2 ' ^ ' 

where T2 and Ti are the postshock and prcshock temperatures. Shock centers identified in 
multiple directional passes were labeled by the maximum M. We followed only those portions 
of shock surfaces with M > 1.5 to avoid confusion from complex flow patterns and shock 
surface topologies associated with very weak shocks. In the actual simulations, the minimum 
gas temperature was set as the temperature of CMBR, that is, Tmm = 2.7(1 + z), since 
photoionization and heating as well as radiative cooling were ignored. However, considering 
that significant reionization would have taken place by z ~ 15 (see, e.g., Haiman & Holder 
2003), in post-processing we set the minimum gas temperature at T^^^ — 10^ K. The shock 
properties described below include that consistency adjustment. 

Figure 1 shows a two-dimensional slice from the 1024^ simulation isolating a typical 
group with X-ray emission weighted temperature, Tj. 1.3 keV, and X-ray luminosity, 
f» 4.2 X 10'^^ h erg s~^. The figure shows the locations of identified shocks along with the 
X-ray emissivity, gas temperature, density, and velocity field distributions. Although the 
X-ray emissivity distribution looks relatively smooth and round, there are complex accreting 
flows around the group including three sheets (dotted line contours in the temperature 
contour map) and one filament (thin solid line contours). A complex topology of shock 
surfaces surrounding the group can be seen in the temperature contour map and through 
the locations of identified shocks. The group as well as the associated sheets and filament are 
bounded by shocks, but there arc several additional shock structures within that represent 
a variety of converging fiow patterns. To illustrate the limitations of conventional spherical 
accretion concepts, we note that this group provides an example of structures where shocks 
can form in the core by organized infall flows accreting from fllamcnts and sheets that 
penetrate deep into the center of potential wells. The Mach number of the surface along 
the portion of the accretion shock centered at {x, y) = (3.5, 4.0)/i~^ Mpc ranges over 2.4 — 4 
with a mean value of 3.2, while that along the portion centered at {x, y) — (5.5, 4.0)/i~^ Mpc 
ranges 5.3 — 8.3 with a mean value of 6.5. Although quite strong shocks of Mach number of 
a few, including these, were often found within 0.5 — lh~^ Mpc from the center of clusters 
and groups, most shocks identifled in our simulations are located outside the cores of clusters 
and groups (see the next section and Figure 5). 
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In order to illustrate how our shock identification scheme works, we plot in Figure 2 

the flow structure along the horizontal path drawn just below y = 4h~^ Mpc in the shock 
location plot of Figure 1. As indicated in the upper right panel of the flgure, two outer shocks 
with high Mach numbers and three inner shocks with lower Mach numbers were identified 
along the path. The existence of these five shocks is also clearly evident in the velocity field 
and temperature contours in Figure 1. Extensive tests showed that our scheme identifies 
shock surfaces reliably with a typical error of a few percent in Mach number. 



3. Properties of Cosmological Shock Waves 

In previous studies by Miniati et al. (2000) and by other authors, cosmological shock 
waves were often organized into three categories; i.e., accretion shocks, merger shocks, and 
internal fiow shocks. However, based on close examination of shock locations and the prop- 
erties of the shocks and their associated flows in our simulations, we suggest instead that 
it is very informative to classify cosmological shocks into two broad populations that can 
be conveniently labeled as external and internal shocks, depending on whether or not the 
associated preshock gas was previously shocked (that is, whether Ti < 10^ or > 10"^ for the 
preshock temperature in practice). This binary classification facilitates the understanding 
of their role in energy dissipation (see the next section). External shocks surround sheets, 
filaments and knots, forming when never-shocked, low density, void gas accretes onto those 
nonlinear structures. Subsequent, internal shocks are distributed within the regions bounded 
by external shocks. They are produced by flow motions accompanying hierarchical structure 
formation inside the bounding shocks. For more reflned questions internal shocks can be fur- 
ther divided into three types: 1) accretion shocks produced by infall from sheets to fllaments 
and knots and from fllaments to knots, 2) merger shocks formed during subclump mergers, 
and 3) flow shocks induced by chaotic supersonic motions inside the nonlinear structures. 

Figure 3 represents a two dimensional slice of a (25/i~^ Mpc)^ volume extracted from 
the 1024^ grid zone simulation. It shows the distributions of external and internal shocks in 
a cluster with X-ray emission-weighted temperature, ^ 3.3 keV, and X-ray luminosity, 
Lx ~ 1.4 X lO^^/i erg s~^, along with the gas density distribution and the velocity fleld 
of the inner parts of the slice. External shocks here define an entire "cluster complex" 
which has dimensions of about (10 x 10 x 20)(/i~^Mpc)^. Numerous internal shocks were 
identified inside the complex. The two-dimensional velocity field demonstrates that infalls 
from several associated filaments and sheets form accretion shocks and also induce chaotic 
flow motions in the medium around the cluster. Figure 4 shows a three-dimensional view 
of shock surfaces around the same cluster complex. Surfaces of external shocks with high 
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Mach numbers, represented with yellow, encompass the complex and filaments associated 
with it. Several sheets with lower Mach numbers, which intersect at the complex, are also 
clearly visible. These figures show clearly that the canonical spherical external accretion 
shock model around a cluster would be far too simple to apply to this sample cluster. Note, 
in particular, that a "spherical" cluster of = 3.3 keV has the first caustic at Rc ~ 3.5h~^ 
Mpc from the center in our ACDM model (Ryu & Kang 1997), around which the external 
shock would be located. That is well inside the external, accretion shock in this simulated 
cluster complex. 

We will now estimate some quantitative measures of shock frequency. To start we com- 
puted the surface area of identified shocks per logarithmic Mach number interval, dS{M, z)/d log M, 
normalized by the volume of the simulation box. This provides an effective inverse comoving 
mean distance between shock surfaces. In this accounting each shock center contributes a 
surface 1.19(AZ)^, which is the mean projected area within a three-dimensional zone for ran- 
dom shock normal orientations. The upper two panels of Figure 5 show dS{M, z)/dlogM 
for external and internal shocks at several epochs in the simulation with 1024^ grid zones. 
Table 1 lists the integrated mean shock separation, 1/S{z), along with the mean values of 
the shock Mach number, the shock speed, the preshock sound speed and the preshock gas 
density. 

Several important points are apparent. 1) The two populations of shocks have distinctive 
distributions of shock surfaces, dS{M, z)/d\og M , implying they arc induced by flows of 
different characteristics. The area distribution of external shock surfaces peaks at M ^ 3 — 5, 
extending up to M ~ 100 or higher. In contrast, the comoving area of internal shocks 
increases to the weakest shocks we identified in our analysis {i.e., M — 1.5). 2) While 
the comoving area of external shock surfaces has not changed much since z < 2, that of 
internal shocks has increased significantly. The former behavior refiects the fact that the 
"bounding" shapes associated with external shocks have evolved towards simpler shapes to 
almost balance the increasing comoving volumes enclosed. On the other hand both the 
volume enclosed and the complexity of internal shocks has increased, so that these shocks 
now include more area. At the present epoch the total area of external shock surfaces is ~ 2 
times of that of identified internal shocks. 3) Although the mean Mach number is higher for 
external shocks, the mean shock speed is actually larger for internal shock. In addition, the 
preshock gas density is significantly higher for internal shocks. These factors enhance their 
dynamical importance (see the next section). 4) We found that most clusters and groups 
with Tj; > 0.1 keV have shocks within 0.5h~^ Mpc from the centers at present. The area 
distribution of these "cluster shocks", shown in the upper right panel of Figure 5, fits best 
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to 



dS{M, z) 
dlogM 



oc exp 




(2) 



with Mch ~ 1 in the range of M < 10. Their mean Mach number is ~ 4. The cluster 

shocks, however, actually account for only a very small fraction of identified shock surfaces. 
We emphasize that the statistics for the cluster shocks would have been affected by the 
finite resolution, Al = 97.7/i"^ kpc, as well as by the exclusion of physical processes like 
radiative cooling and feedback from galaxies and stars that infiuence conditions inside cluster 
cores. Still, it is significant that, compared with the distribution of binary merger shocks 
studied in Gabici & Blasi (2003), shocks with higher Mach numbers are more common in the 
environments of clusters and groups in our simulation. This difference is because our cluster 
shock population includes accretion shocks created by infall from filaments and sheets to 
knots and fiow shocks generated by chaotic supersonic motions, as well as those induced by 
hierarchical merging. Especially, internal accretion shocks are strong with a typical Mach 
number of several (see §2.2). 

The lower two panels of Figure 5 show dS (M, 0)/d log M for external and internal shocks 
at z = in the simulations with different numerical resolutions of 1024"^ — 128^ grid zones. 
Table 2 lists the integrated surface area, S{0), along with the mass of the gas that went 
through shocks at least once, ^^^^(O), and the total gas thermal energy inside the compu- 
tational box, Sth{0), at the present epoch in the simulations with 1024^ — 64^ grid zones. 
The shocked gas mass, J^sg was estimated by summing the mass of gas with T > lO'' K. 
We expect that finite resolution would affect the statistics of shocks in two different ways: 
1) some shocks, especially weak ones with low Mach numbers, may not have been "cap- 
tured" in the hydrodynamic part of simulations, and 2) they may not have been "identified" 
by our post-processing shock identification scheme, especially in the regions with complex 
three-dimensional fiow structures. The former causes quantities like Aigg and £th to be un- 
derestimated. On the other hand, the latter reduces the estimated value of S{z) and so the 
estimation of shock dissipation (see the next section). We see that Aigg and Sth converge 
rather quickly and their converged values should be within ~ 10% or so of those from the 
simulation with 1024^. This indicates that the abihty of our code to capture shocks is not 
the major limitation. The convergence of the mean shock separation, 1/5" and Sext/Smt, is 
not as obvious from the table. The difficulties in reaching final convergence in these quan- 
tities is evident in Figure 5. Detection of strong shocks with M > 10, especially strong 
external shocks, is relatively robust in the simulation data with 256^ or more grid zones. 
However, the area spanned by weak shocks continues to increase with resolution, reflecting 
the increased flow complexity captured inside nonlinear structures. The quantities of l/S" 
and Sext/ Sint are, however, also plotted in Figure 8, where it is evident that convergence is 
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underway between the two highest resolution simulations. It is likely that l/S" ~ 4/i~^ Mpc, 
and Sext/Sint ~ 2. It is, then, interesting to note that since the volume of the nonlinear 
structures where internal shocks are found is ~ 1/10 of the entire computed volume, the 
mean distance between internal shock surfaces is ~ Ih^^ Mpc within the nonlinear struc- 
tures, which is comparable to the scale of the nonlinear structure involved. This conclusion, 
of course, includes the caveats that only shocks with M > 1.5 have been counted and that we 
have omitted from the simulations physics likely to influence shock formation inside cluster 
cores. Accepting the above limitations, we argue in the next section from estimations of the 
gas thermal energy dissipated at shocks and the mass passed through external shocks that 
the identification of dynamically important shocks should be complete within an error of 
order 10% in the simulation data with 1024^ grid zones. 



As the first step to quantitative estimates of the gas thermal energy and CR energy 
from dissipation at cosmological shock waves, we defined the following fluxes of mass and 
energies at each "shock center": 1) the gas mass flux incident on shocks, fms {— PiVsh)', 2) 
the incident kinetic energy flux, (= (l/2)piv^^); 3) the thermal energy flux generated at 
shocks, fth (defined below); and 4) the CR energy extracted, i.e., "nonthermal dissipation", 
at shocks, fcR (defined below). Here subscripts 1 and 2 stand for preshock and postshock 
conditions, respectively. 

We define thermal energy fiux generated at shocks as 



We point out that the second term inside the brackets subtracts the effect of adiabatic 
compression occurred at a shock, not simply the thermal energy fiux entering the shock; 
namely, eth,iUi. The ratio fth/ f(t> = S{M) then defines the efficiency of shock thermalization, 
which is a function of Mach number only, and can be determined from the Rankine-Hugoniot 
jump conditions. The left panel of Figure 6 shows S{M). As expected, S{M) increases with 
Mach number, asymptotically approaching 5{M) — > 0.56 for M > 10. 

We express the CR energy extraction rate as fcR — r}{M)f^, where r]{M) measures 
the efficiency of diffusive shock acceleration for a given Mach number. To estimate this 
efficiency we have used results of numerical simulations of the nonlinear evolution of CR 
modified quasi-parallel plane shocks (e.g., Kang et al. 2002; Kang & Jones 2002). The 
simulations utilize a plasma-physics-based model for the leakage of thermal ions into the CR 



4. Energy Dissipation at Shocks 




(3) 
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population, which depends on a single parameter measuring the enhancement of nonhnear 
scattering waves across the shock transition. That parameter is reasonably well limited by 
the theory and by plasma simulations (see, e.g., Malkov 1998), and, in the DSA simulation 
we used, typically limits the fraction of protons injected into the CR population at shocks 
to be of order ~ 10^'^. The CR acceleration efficiency characteristics of our DSA model 
are broadly consistent with other widely used theoretical and numerical studies of nonlinear 
CR shocks (e.g., Malkov 1999; Berezhko et al. 1995; Ellison et al. 1996). In order to apply 
properties of our DSA model to cosmological shock waves, we calculated proton acceleration 
and accompanying CR- modified fiow evolution for shocks with Vgh — 1500 — 3000 km 
propagating into media of Ti = 10^ — 10^ K, assuming Bohm-type diffusion for the CRs (for 
details see Kang et al. 2002). As the shock structures evolved we determined the ratio of 
the total CR energy extraction during the evolution of a shock to the kinetic energy passed 
through the shock according to 



The values of $cij(^) quickly reached approximate "time-asymptotic" values, after which the 
shock structures evolved approximately "self-similarly." These asymptotic values of ^cR{t) 
were taken as our estimates for rj{M). 

The left panel of Figure 6 shows the resulting CR acceleration efficiency, ri{M), along 
with the gas thermalization efficiency 5{M). Both 5{M) and 7;(M) increase with Mach 
number. In strong, high Mach number shocks, r]{M) approaches the asymptotic value of 
ri{M) — >~ 0.53. We see that 5{M) is somewhat larger than ri{M) over all Mach numbers. 
On the other hand, we note that 5{M) was computed for purely gasdynamical shocks; that 
is, without accounting for energy removal into CRs and enhanced adiabatic compression 
within the CR shock precursor. If nonlinear dynamical feedback of CRs were included self- 
consistently in estimating 5(M), the resulting 5{M) would be somewhat smaller than that 
shown in Figure 6. 

From the above gas mass and kinetic energy fluxes at each shock center, we calculated 
the associated fluxes through surfaces of shocks with Mach number between logM and 
logM-|-(i(logM) at different redshifts; i.e., dFj^g^M, z) and dF^{M,z), respectively. We also 
calculated the similarly defined fluxes of gas thermal and CR energies dissipated at shocks as 
dFth{M, z) = dF^{M, z)xS{M) and dFcR^M, z) = dF^{M, z)xr]{M), respectively. The right 
panel of Figure 6 illustrates dF^{M,z)/dlogM per unit comoving volume for external and 
internal shocks at different redshifts in the simulation with 1024^ grid zones. A noticeable 
point is that the kinetic energy flux, and also the gas mass flux through external shocks were 
larger in the past. This is because the preshock gas density was larger in the past. However, 





the kinetic energy flux through internal shocks has been more or less constant since z ~ 1.5, 
and was smaller before that. 



To provide measures of the roles of shocks over time we integrated from z — 2 to z — 
the gas mass and kinetic energy that passed through shock surfaces and the gas thermal and 
CR energies dissipated at shock surfaces as 

dYi{M) _ 1 ['=UF,{M,z\^^ 



dlogM MJz=2 d\ogM 

where the subscript i = ms , (f), th, and CR stands for the four fluxes deflned above. The 
quantities were normalized either to the shocked gas mass, Mm.s = M.sg{z = 0), for mass or 
to the total gas thermal energy inside the computational box, A/i = Sth{z = 0), for energies 
(see Table 2). We also summed these time- integrated measures to calculate the associated 
global shock-processed quantities. 



Yi{> M) 



I'M 


'dYi{My 


1 oo 


_d\ogM_ 



dlogM. (6) 



Yms and dYms{M)/d\ogM not only measure the total mass that passed through shocks in 
the z = 2 to z = interval, but also, by way of their normalization, measure the mean 
number of times that gas has been subjected to shock dissipation. Meanwhile, Y^ and its 
derivative measure the total kinetic energy that has been subject to shock dissipation. On 
the other hand, Yth and Ycr compare the total thermal and CR energies that have resulted 
from shock dissipation. 

Figure 7 shows dYi{M)/dlogM and Yi{> M) for external and internal shocks in the 
simulation with 1024^ grid zones. Again, several important points are apparent. 1) The plots 
for dYms{M)/ dlogM and Yjns{> M) indicate that more mass has passed through internal 
shocks than external shocks. With Yms{> 1-5) ^ 2.2 for internal shocks, the mass inside the 
nonlinear structures of sheets, fllaments and knots has been shocked, on average, twice or so 
by internal shocks from 2; = 2 to 0. But we note that dYjns{M) / dlog M for internal shocks 
increases to M = 1.5, the lowest Mach number we kept, and weak shock identiflcation is 
not fully converged yet with this numerical resolution. So the above value of 2.2 should be 
regarded as the lower limit. 2) On the other hand, since gas enters the nonlinear structures 
by passing through external shocks, l^m.s(> 1-5) for external shocks {i.e., the mass passed 
through external shocks from 2; = 2 to 0) should match the increase in the shocked gas mass, 
Msg: from 2; = 2 to 0. In the simulation with 1024^ grid zones, we get Yjns{> 1-5) = 0.35 for 
external shocks, while 0.42 = (0.73 - 0.42)/0.73 is expected with Msg{z = 2) = O.A2Mgas 
and M.sg{z = 0) = 0.73M.gas (where M.gas is the total gas mass in the computational 
box). Although a fraction of the hot gas may have been heated above 10^ K by adiabatic 
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compression in void regions or by weak shocks with M < 1.5, some of this discrepancy could 
be due to the under-counting of external shocks. So we estimate an error of ~ 10% or so in 
the completeness of the identification of external shocks. 3) The lower two panels of Figure 7 
show that internal shocks play a more important role than external shocks in the dissipating 
energy associated with structure formation. Specifically, internal shocks with 2 < M < 4 
account for ~ 1/2 of dissipation. While the thermal energy generation peaks for shocks in 
the range 2 < M < 3, the CR energy extraction peaks in the range 3 < M < 4. Miniati 
(2002) identified shocks with 4 < M < 5 as the most important in gas thermahzation. 
Our result indicates a somewhat lower Mach number range, because we found more weaker 
shocks in our simulation data. This difference is probably due to our greater resolution (see, 
e.g.. Figure 5), as well as due to improved shock capturing in the present code and a more 
sophisticated shock identification scheme. 4) The amount of thermal energy generated at 
shocks from 2; = 2 to was computed to be ~ 1.5 times the current gas thermal energy 
in the simulation with 1024^ grid zones. Since no other heating or cooling physics was 
included, there are two likely contributors to the difference. First, we attribute some of 
the difference to adiabatic expansion of nonlinear structures after formation. Expansion of 
~ 8.5% along each dimension would be enough to reduce the gas thermal energy by a factor 
of 1.5. In addition, some of the discrepancy could be the result of multiple counting of 
internal shocks in the regions with complex flow structures. So again, we put an error of 
a few 10%, at most, in the completeness of the identification of shocks. 5) With the DSA 
model we adopted (Kang et al. 2002; Kang & Jones 2002) the ratio of the CR to gas thermal 
energies dissipated at cosmological shock waves with the Mach number greater than 1.5 is 
Ycr{^ l-5)/5^t/i(> 1.5) ~ 1/2. 6) There are many shocks with M < 2, but they are not 
important in energetics. Hence, although we ignored the shocks with M < 1.5 in this work, 
that should not significantly impact quantitative energetics results presented here. 

Figure 8 shows Kj(> 1.5) through all shock surfaces in the simulations with different 
resolution of 1024'^ — 64'^ grid zones. Note that Fj's were normalized with Msg and 8sg for a 
given resolution (Table 2). Yet, for instance, F^, the integrated kinetic energy flux through 
shock surfaces, decreases in the plot by 60% from 1024^ to 512^ and by 95% from 512^ to 
256^. On the other hand, Egg-i the total gas thermal energy inside the computational box, 
decreases only by 8% from 1024^ to 512'^ and 13% from 512^ to 256^. This reinforces the 
statements in the previous section that the error in quantitative assessments of this paper 
comes mostly from shock identification in the post-processing analysis rather than shock 
capturing in the code. This also shows that our estimation of an error of order 10% in the 
data of the 1024^ simulation is consistent with resolution convergence. 
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5. Summciry 

We identified and studied shock waves with Mach number M > 1.5 in a set of cosmo- 
logical N-body/hydrodynamic simulations for a ACDM universe in a cubic box of comoving 
size 100/i^^ Mpc. To facihtate the analysis of their properties, the cosmological shock waves 
were classified as external and internal shocks, depending on whether or not the preshock gas 
was previously shocked. External shocks form around outermost surfaces that encompass 
nonlinear structures, so they are by nature accretion shocks that decelerate the previously 
unshocked intergalactic gas infalling toward sheets, filaments and knots. Internal shocks are 
produced within those nonlinear structures by accretion fiows of previously shocked gas from 
sheets to filaments and knots and from filaments to knots, by merging of subclumps, or by 
chaotic flow motions induced in the course of hierarchical clustering. 

For all shocks of Af > 1.5 identified in the simulation of highest resolution, the mean 
distance between shock surfaces, the inverse indicator of shock occurrence, is ~ 4h-^ Mpc at 
present. Further, external shocks arc more extensive, with their surface area ~ 2 times larger 
than that of identified internal shocks at present. With the volume of nonlinear structures 
which is ~ 1/10 of the total volume, the mean distance between internal shock surfaces 
is ~ lh~^ Mpc within nonhnear structures at present. Although external shocks typically 
have higher Mach numbers, internal shocks have higher shock speed and higher preshock gas 
density. As a result, internal shocks are responsible for ~ 95% of gas thermahzation and for 
~ 90% of CR acceleration at shocks, and they process about 6 times more gas mass through 
shock surfaces than external shocks do from z — 2to z = 0. Internal shocks with 2 ^ Af < 4 
are especially important in energy dissipation, contributing ~ 1/2 of the total. By adopting 
a model of CR proton acceleration based on nonlinear diffusive shock simulations (Kang et 
al. 2002), our study predicts that the ratio of the CR to gas thermal energies dissipated at 
all cosmological shocks through the history of the universe could be ~ 1/2. Due to long CR 
proton trapping times and energy loss lifetimes, they should fill the volumes inside filaments 
and sheets as well as in clusters and groups. Short electron lifetimes, however, lead that 
population to depend on other factors (Miniati et al. 2001a). The existence of substantial 
CR populations could have affected the evolution and the dynamical status of the large scale 
structure of the universe (see, e.g., Enfilin et al. 1997). 

The examination of results from simulations of different resolutions showed that the 
convergence with resolution in shock identification is slower than in shock capture within 

the simulation itself. Consistency checks and resolution convergence analysis lead to error 
estimates of order of 10% in our quantitative estimates of the accounting of and energy 
dissipation in cosmological shock waves. 

Strong, external shocks are energetically less important than moderate strength internal 
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shocks. However, it is important to keep in mind that the large curvature radii and long life 
times of these shocks make them viable candidates to accelerate CRs to ultra high energies 
of several xlO^^ eV (Norman et al. 1995; Kang et al. 1996). 

The Mach number distribution and the amount of energy dissipation at cosmological 
shocks have significant implications for several cosmological observations such as radio and 
7-ray emissions as well as X-ray emission from the ICM and contribution to the cosmic 
7-ray background from CRs accelerated at these shocks (see, e.g., Loeb & Waxmann 2000; 
Miniati et al. 200 lb, a; Miniati 2002). These important issues, along with the observational 
manifestation of cosmological shock waves, will be considered in future studies. 
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Table 1. Mean flow quantities of external/internal sliocks at several different epochs 



z 


IjS ^ 




(M)ext 




{Vsh)ext 


" {Vsh)mt 


a /„ \ a 
\"--s/ext 


(Cs)int 


^ {Psh)ext 


^ (Ps/i)iE 





4.4 


2.1 


8.0 


3.2 


123 


226 


15.3 


82 


1.05 


6.78 


0.2 


4.4 


2.3 


8.1 


3.3 


123 


230 


15.3 


83 


1.12 


7.15 


0.5 


4.5 


2.8 


8.0 


3.3 


122 


231 


15.3 


83 


1.25 


7.86 


1 


5.0 


3.7 


7.5 


3.4 


114 


214 


15.3 


76 


1.48 


8.87 


1.5 


5.7 


5.0 


7.0 


3.4 


107 


196 


15.3 


69 


1.79 


10.3 


2 


6.8 


6.6 


6.5 


3.4 


100 


177 


15.3 


62 


2.14 


10.9 



^Lengths in units of (1 + ^Mpc, speeds in km s ^, and density compared to tlie 

mean comoving density of gas {Pgas){z), respectively. 



Table 2. Shock Associated Quantities Measured with Different Resolutions 



resolution 


1/S ^ 


'S'ext/ 'S'int 


Msg ' 


c a 


1024^ 


4.4 


2.1 


0.73 


1.42 


512-'^ 


7.3 


3.2 


0.69 


1.30 


256^ 


14 


5.4 


0.61 


1.13 


128=^ 


39 


9.6 


0.45 


0.87 


643 


260 


28 


0.19 


0.40 



^Lengths in units of h ^Mpc, mass compared to total gas mass inside the computational 
box Mgas, and energy in units lO^h"^ ergs, respectively. 
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Fig. 1. — X-ray emissivity contours (upper left panel), velocity field superimposed on density 
contours (upper right), gas temperature contours (lower left) and shock locations (lower 
right) in a two-dimensional slice of {9.76h~^ Mpc)^ around a group of X-ray emission weighted 
temperature ^ 1.3 kcV at 2; = 0. Contours of gas density Pgas/{Pgas) > 1 ai'c shown. In 
the temperature contours, heavy solid lines are used for logT > 6.8, light solid lines are for 
5.5 < logT < 6.8, and dotted lines are for 4 < logT < 5.5. Several sheets and a filament 
were identified in the temperature contours. 
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Fig. 2. — Flow structure along the line path drawn in Figure 1, showing gas density, Mach 
number of identified shocks, gas temperature and V • v. 
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Fig. 3. — Two-dimensional slice of (25/i~^Mpc)^ around a complex including a cluster of 
X-ray emission weighted temperature ~ 3.3 at z = 0, showing gas density, internal and 
external shock distributions, and velocity field. For clarity, the velocity field is shown in a 
zoomed region of (12.5/i~^ Mpc)^ centered at the cluster. 
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Fig. 4. — Three-dimensional shock surfaces in a volume of (25/i ^Mpc)^ around the same 
complex as in Figure 3. The color bar shows the values of Mach numbers of shock surfaces. 
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Fig. 5. — Top two panels: Inverse of the mean comoving distance between shock surfaces 
with Mach number between logM and logM + d(\ogM) at different redshifts, dS{M,z), 
for external shocks and internal shocks in the simulation with 1024^ grid zones. The curve 
labeled by "clusters" shows the quantity for shocks inside clusters and groups at present, 
which was multiplied with 5 for clarity (see the text for details). Bottom two panels: The 
same quantity at 2; = in the simulations with different resolutions of 1024^ — 128^. 
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Fig. 6. — Left panel: Gas thermalization efficiency, S{M), and CR acceleration efficiency, 
r]{M), at shocks as a function of Mach number. Dots for 77 (M) are the values estimated from 
numerical simulations based on a DSA model and solid line is the fit. Right panel: Kinetic 
energy flux per unit comoving volume through surfaces of external and internal shocks with 
Mach number between logM and logM + d{\ogM) at different redshifts, dF^{M, z), in the 
simulation with 1024^ grid zones. The line types are same as those in the upper panels of 
Figure 5 
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Fig. 7. — Upper left panel: Mass, dYms{M), processed through surfaces of external and 
internal shocks with Mach number between logM and logM + rf(logM), from 2; = 2 to 
z = 0. Upper Right panel: Mass, Ymass{> M), processed through surfaces of external and 
internal shocks with Mach number greater than M, from z = 2 to z = 0. Lower left panel: 
Kinetic energy, dY^{M), thermal energy, dYth{M), and CR energy, dYcR{M), processed 
through surfaces of external and internal shocks with Mach number between logM and 
logM + (i(logM), from 2; = 2 to 2; = 0. Lower Right panel: Same energies, Y^{> M), 
Yth{> M), Ycr{> M), processed through surfaces of external and internal shocks with Mach 
number greater than M, from 2; = 2 to 2; = 0. The mass and energies are normalized to the 
the shocked gas mass, Aigg, and the total gas thermal energy, Sth, inside simulation box at 
z — 0, respectively. All plots are from the simulation data with 1024^ grid zones. 
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Fig. 8. — Left panel: Mean shock separation ed, z — 0, 1/S, and its ratio for two shock 
populations, Sext/Sint, for identified shocks with M > 1.5 in the simulations of different 
resolutions with number of grid points along one-dimension, — 64, 128, 256, 512, and 
1024. Right panel: Integrated mass, Yms, and energies, Y^t,, Yqr, processed through all 
identified shocks with M > 1.5 in the simulations of different resolutions. The mass and 
energies arc normalized to the shocked gas mass, A4sg- and the total gas thermal energy, Sth, 
of given resolution, inside simulation box at z = 0, respectively. 
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